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ABSTRACT 



Modelling the formation of super-km- sized planetesimals by gravitational collapse of regions overdense in small particles requires 
numerical algorithms capable of handling simultaneously hydrodynamics, particle dynamics and particle collisions. While the initial 
phases of radial contraction are dictated by drag forces and gravity, particle collisions become gradually more significant as filaments 
contract beyond Roche density. Here we present a new numerical algorithm for treating momentum and energy exchange in colli- 
sions between numerical superparticles representing a high number of physical particles. We adopt a Monte Carlo approach where 
superparticle pairs in a grid cell collide statistically on the physical collision time-scale. Collisions occur by enlarging particles until 
they touch and solving for the collision outcome, accounting for energy dissipation in inelastic collisions. We demonstrate that super- 
particle collisions can be consistently implemented at a modest computational cost. In protoplanetary disc turbulence driven by the 
streaming instability, we argue that the relative Keplerian shear velocity should be subtracted during the collision calculation. If it is 
not subtracted, density inhomogeneities are too rapidly diff'used away, as bloated particles exaggerate collision speeds. Local particle 
densities reach several thousand times the mid-plane gas density. We find efficient formation of gravitationally bound clumps, with 
a range of masses corresponding to contracted radii from 100 to 400 km when applied to the asteroid belt and 150 to 730 km when 
applied to the Kuiper belt, extrapolated using a constant self-gravity parameter. The smaller planetesimals are not observed at low 
resolution, but the masses of the largest planetesimals are relatively independent of resolution and treatment of collisions. 
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1. Introduction 



methods: numerical - minor planets, asteroids: general - planets and satellites: formation - protoplan- 



The formation of super-km- sized planetesimals is an im- 
portant step towards terres trial plane t s and the solid cores 
of gas and ice giants (e.g. ISafronovi 1 19691: iGoldreich et all 
120041: IChiang & YoudinL |201Q|) . The asteroid and Kuiper belts 



of the solar system, as well as the extrasolar debris discs, 
are believed to be left-over populations of planetesimals that 
did not grow to planets. Comparing models and simulations 
of planetesimal formation to observations of such planetesi- 
mal belts constrains our theoretical picture of the planetes- 
imal formation stage, and at the same time it gives insight 
into the physical processes that shaped the architectures of 
these systems ( [Morbidelli et al., 2009; Weidenschilling, 2010; 

m 



Nesvornv etalir2 QlQ: She ppard & Truiilldl201Ql:IKrivovLI20ia: 
Kenvon & Bromley. 201 Qi) . 



Planetesimal formation takes place in a complex environ- 
ment of turbulent gas interacting via drag forces with parti- 
cles of many sizes. The streaming instability thrives in the 
systematic relative motion of gas and particles and leads to 
spontaneous clumping of particles ( Youdin & Goodmanl l2005l: 



^ pin^ pa 

iJohansen & YoudinL l2QQ7l: iBai & StoneLl2QlQbl). seeding a g rav- 
itational collapse into bound clumps fJohansen et aP. 120091) and 



further to solid planetesimals (jNesvorny etaLj, 12010D. While 
the latest years have seen major progress in numerical mod- 
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elli ng of drag forc e inte r action between particles and gas 
(Youdin &Johanseni 120071: iBalsara et al.L l2009l: iMiniatiL I201QI: 
iBai & Stone. '201 0al) as well a s the self-g r avity of the parti- 
cle layer (Johansen et al.L l2007l: iRein et all l2010l) . good algo- 
rithms for treating simultaneously hydrodynamics, gravitational 
dynamics and particle collisions are still missing. 

There are two main approaches in astrophysics to treat- 
ing particle collisions in numerical simulations. Modelling a 
set of physical particles with collision tracking allows simu- 
lation of particle aggregation in close concordance with the 
nature of real physical collisions. This method has success- 
fully _been_a22lied to inod parti c le rings of Saturn 
(Wis dom &Tremainel Il988l: ISalol 1 19911: iKarialainen & Said 
>2QQ4l) and to n iodel collisions between i ndividual dust grains 
and aggregates (iDominik & Nuboldil2002l) . The drawback of the 
physical-particle approach is that the size of the system is limited 
by the number of numerical particles that can be afforded in the 
simulation. The formation of a Ceres-mass planetesimal from 
10-cm- sized rocks would e.g. require tracking of (9(10^^) par- 
ticles, orders of magnitude beyond what current computational 
resources allow. 

Algorithms involving inflated particles group collections of 
physical particles into much larger numerical particles under 
conservation of total mass M and mean free path X. Decreasing 
the particle number A/^ to a number that can be handled in a com- 
puter simulation, while maintaining = (N/V)cr by artificially 
increasing the collisional cross section cr, yields the correct col- 
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lision frequency in systems that are much larger than what can 
be resolved with the physical particle approach. The inflated par- 
ticle approach was used recently by, Lithwick & C hiang ( 20071 
iMichikoshi et alJ (l2007h . iNesvornv et al.lT201Ql) . and lRein et ^ 
(l201Qh . with different methods for tracking the actual colli- 
sion, but the concept of bloated particles has deeper roots (e.g. 
lKokubo&IdaLll996h . 

In this paper we put forward a new algorithm to model col- 
lisions between numerical sup erp articles. Superparticles are de- 
signed to represent swarms of physical particles. The aerody- 
namical properties of the superparticle (e.g. the friction time) is 
still that of a single physical particle. Superparticles are widely 
used to model the solid particle component in computer sim- 
ulatio ns of coupled g as and particle motion in protoplanetary 
discs (iJohansen & Yo udin, 2007; Bai & Stone, 2010b). Since su- 
perparticles can be considered to represent swa rms of smaller 
particl es, direct collision tracking is not possible, [johansen et al.l 
(l2007l) modelled superparticle collisions by damping the random 
motion of particles inside a grid cell on the collisional time- 
scale. They showed that inelastic collisions, where part of the 
kinetic energy is converted to heat and deformation during the 
collisions, is beneficial for the gravitational collapse and allows 
the formation of planetesimals in protoplanetary discs of lower 
mass, compared to simulations without damping. However, the 
simplified collision scheme of Johansen et al. ( 2007) is insuffi- 
cient in capturing the pairwise momentum exchange and energy 
dissipation. 

We develop here a statistical approach to model the full mo- 
mentum exchange and energy dissipation in collisions between 
superparticles. The Monte Carl o scheme is inspired by the col- 
lision algorithms pr e sented by iLithwick & Chiang! (l2007b and 
IZsom & Dullemondl (l2008l) . The essence of our algorithm is to 
determine the collision time- scale between all superparticle pairs 
within a grid cell. Two superparticles collide as if they were 
physical particles touching each other, if a random number cho- 
sen uniformly between zero and one is smaller than the ratio of 
the simulation time- step to the collision time- scale. 

Collisions can be followed together with hydrodynamics at 
a moderate computational cost depending only on the number of 
particles per grid cell. We compare the statistical properties of 
the particle density in 3-D hydrodynamical simulations with and 
without collisions. Including the self-gravity of the particles, we 
find formation of bound clumps, with masses comparable to that 
of the 500-km-radius dwarf planet Ceres when applied to the as- 
teroid belt, relatively independently of numerical resolution and 
treatment of collisions. The scale-free nature of our simulations 
allows application of the results to the Kuiper belt as well, with 
contracted planetesimal radii approximately 80% higher than in 
the asteroid belt. 

The paper is organised as follows. In Sect. |2] we describe the 
new superparticle collision algorithm. The algorithm is tested 
against known test problems and conservation properties of the 
shearing box in Sect. O In Sect. IH we analyse statistical prop- 
erties of the particle density achieved in simulations of gas and 
particle turbulence driven by the streaming instability. We con- 
tinue to include self-gravity in the simulations and analyse the 
planetesimal masses obtained under various assumptions about 
collisions in Sect. [51 We summarise and discuss our results in 
Sect. [6j The appendices A-C contain further descriptions of the 
collision algorithm. 



2. Superparticle collision algorithm 

We will use the notation that a superparticle represents a swarm 
of physical particles with number density h and volume 6V. 
Since we are interested in coupling superparticle collisions to 
grid hydrodynamics, the volume is taken to be that of a grid cell, 
dV = dxxdy X Sz. The physical particles in the swarm have in- 
dividual mass, physical radius, material density, and collisional 
cross section m, R, p, and cr. We assume that all swarms are sim- 
ilar, both in internal particle number and in the physical mass of 
the constituent particles. 

To track a collision we calculate the mean free path A for a 
test particle interacting with the swarm of particles represented 
by a single superparticle. 



1 

^= — . 

ncr 



(1) 



Superparticles in the same grid cell are considered as potential 
colliders. For each collision pair the collision time- scale is cal- 
culated from 



A 

Sv' 



(2) 



where Sv is the relative speed between particles / and j. The sim- 
ulation time-step 6t, set by hydrodynamics and drag forces, is 
then used to calculate the probability that those two particles 
collide in this time-step. 



Tc 



(3) 



Two colliding swarms have their velocity vectors changed 
instantaneously. The collision outcome is found by considering 
two virtual spherical particles whose surfaces touch, with par- 
ticle centres at the locations of the superparticles, and solving 
for momentum conservation and inelastic energy dissipation (or 
energy conservation, in case of elastic collisions). We define the 
velocity vectors relative to the mean velocity field v = (Vj-\-Vk)/2, 



= Vk-V = -V . 



(4) 
(5) 



Here vj and Vk are the velocity vectors of the two particle^ The 
normal vector e_i connecting the centres of the particles at the 
time of collision is calculated as 



Xj AT^ 

\Xj-Xk\ 



(6) 



The parallel vector is perpendicular to in the same plane as 
the relative velocity vector. The relative velocity vectors are now 
decomposed on the two directions 



(7) 
(8) 



with Gk = -aj and bk = -bj. In the collision we maintain b, 
while we reflect a according to 



a 



-ea . 



(9) 



^ We show in Sect. 13.2.11 that the Keplerian shear should be sub- 
tracted from the velocity vectors when determining both the collision 
time- scale and the collision outcome, in the limit of particles that are 
much smaller than a grid cell. 
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Here e e [0, 1] is the coefficient of restitution, parameterising the 
degree of energy dissipation during the collision. Inelastic col- 
lisions can play an important role in dissipating kinetic energy 
and facilitating the gravitational collapse phase. In general the 
coefficient of restitution depends on material parameters, impact 
speed and ambient temperature. Water ice particles have been 
measured to have a high coefficient of restitution e ^ 0.9 for im- 
pact s peeds below ^ 2 m/s (quasi-elastic regime of iHiga et al.L 
Il996h . Above this critical speed the measured coefficient of resti- 
tution rapidly drops towards zero. More recent microgravity and 
drop tower experiments find a coefficient of restitution between 
0.06 and 0. 84 in low-velocity collisi ons between 1.5-cm-sized 
icy pebbles (iHeiBelmann et al.Ll2010l) . In this paper we consider 
for the sake of simplicity the coefficient of restitution to be a 
constant that is independent of the relative speed. 

The collision time- scale has a simple relation to the friction 
time- scale when particles are small and drag forces are in the 
Epstein regime. We show in Appendix lAl how the collision time- 
scale can be easily calculated from the friction time- scale, useful 
e.g. for simulations of gas and particles in protoplanetary discs. 

Consider now a grid cell containing N superparticles. For 
particle / the collision probability for a representative particl^l 
from superparticle / to collide with the particle swarms j = i -\- I 
to j = N is calculated. The collision occurs if a random num- 
ber, drawn for each collision partner, is smaller than P from 
Eq. ©. The collision instantaneously changes the velocity vec- 
tors of both particles / and j. This way the correct collision fre- 
quency is obtained for both particles, even though the algorithm 
only considers the possible collision / with j, but not j with /. In 
Appendix|B]we describe how to consistently limit the number of 
collision partners, and thus save computation time, in grid cells 
which contain many (» 100) particles. 

There are several advantages to using such a probabilistic 
swarm approach to particle collisions. We mention here a few: (i) 
it is fast because we do not have to track when particles touch or 
overlap within the grid cells, (ii) it allows us to freely choose the 
relative speed that enters the collision frequency, useful e.g. for 
subtracting off' the Keplerian shear (see Sect. 13.2. it , and (iii) the 
algorithm is easily generalisable to also include a probabilistic 
approach to particle coagulation and shattering. 

In Fig. [T] we show the collision path length of test particles 
injected into a medium with 10 superparticles per grid cell and 
a mean free path of A = 0.1. Collisions are tracked through the 
Monte Carlo method described above. The collision algorithm 
makes some particles collide after a short ffight path and others 
after a longer. The distribution plotted in Fig. [T] follows closely 
the expectation N = Nq exp(-£/A). The Monte Carlo approach 
to collisions is very similar to the physical particle approach in 
the distribution of free ffight paths. 

The main technical diff'erence between using inffated parti- 
cles (see introduction) and our newly developed collision algo- 
rithm for superparticles is that inffated particles always collide 
when they overlap physically (the particle size can be associated 
with the grid cell size), while superparticles sharing the same 
grid cell collide with a certain probability which guarantees that 
collisions occur on the average after a collisional time- scale. 
Another diff'erence is that superparticles which do not approach 
must still be allowed to collide, as otherwise the mean free path 
will be too long. Non-approaching particles are collided by ffip- 
ping the relative velocity vector before collision and reffipping 



^ IZsom & DullemonS (l2008h define a representative particle from a 
swarm as a test particle (a random particle from the swarm) used to 
probe the collision time- scale with another swarm. 




Free path/?i 

Fig. 1. Cumulative free path for 1000 superparticles released into 
medium with mean-free-path of A = 0.1. The distribution func- 
tion follows the analytical expectation N = Nq exp(-£/A) very 
closely. Our Monte Carlo algorithm for superparticle collisions 
gives a free path in good agreement with the real physical system 
consisting of many more particles. 

afterwards. The main issue with approaching collisions is that 
collisions occur in fixed grid cells which are not centred on the 
superparticle in question, and thus a superparticle at the edge of 
a grid cell will have too few collision partners if only approach- 
ing collisions are allowed. We show in Appendix ICl how the su- 
perparticle approach transforms smoothly to the inffated particle 
approach when the number of superparticles is reduced. 

The Monte Carlo collision scheme presented here could 
equally well be formulated in terms of inffated particles, by con- 
structing inffated particles smaller than a grid cell. Solving sta- 
tistically for the collision outcome of these "sub-grid" particles 
is mathematically equivalent to the interpretation, chosen for this 
paper, of the numerical particles as swarms. 

3. Validation of algorithm 

We have implemented the Monte Carlo superparticle colli- 
sion scheme described in Sect. [2] into the open source code 
Pencil Cod^ The Pencil Code evolves gas on a fixed grid 
and has fully parallelised modules for a n additional solid com- 
ponent represented by s uperparticles (I Johansen et al.L 120071: 
lYoudin & Johanseni l2007l) . We first validate the collision algo- 
rithm in the limit of inffated particles (i.e. where two particles 
occupying the same grid cell always collide and only approach- 
ing collisions are consider ed), to corn pare our results directly 
to those of IQt hwick & C hiangI (l2007h . The 2-D algorithm of 
iLithwick & ChiangI (120071) has a probabilistic approach to deter- 
mine whether two particles are in the same vertical zone when 
they overlap in the plane. Their algorithm can thus be seen as 
a hybrid of the inffated particle approach and a Monte Carlo 
scheme. 

We set up a test prob lem similar to the one presented in 
iLithwick & ChiangI (l2007b . We define a 2-D simulation box cov- 

^ The code, including the developments de- 
scribed in this paper, can be freely downloaded at 
,http : //code . google . com/p/pencil-code/i 
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ering the spatial interval [-2, +2] x [-2, +2] with 4000 grid cells 
in both the x and y direction. 10"^ particles are placed randomly 
in a ring of full width 0.08 centred at the radial distance r = 1. 
A central gravity source, of strength GM = 1, is placed in the 
centre of the coordinate frame. 

We integrate the particle orbits, including collisions, for 10"^ 
revolutions of the ring centre. In order to compare directly with 
iLithwick & Chiand 0007.) we use their 2-D approximation. The 
particle number density can be approximated as n ~ where 
E is the column (number) density and H is the scale height of 
the particle disc. The random particle motion u can be written as 
u ~ HQ. This yields a collision time 



J2D) 



EcrQ 



T 



(10) 



where r = Zcr is the vertical optical depth of the disc and 
Torb = 2nlQ is the orbital time-scale. While the collision time- 
scale in general depends on the random particle motion, this 
dependence vanishes in the 2-D Keplerian disc approximation 
- faster random motion cancels with increased particle scale- 
height in the collision time expression. 

Requiring that orbits are maintained for 10"^ orbital time- 
scales, we set the time- step of the Pencil Code to = 0.01i2"\ 
covering each orbit Torb = 27ilQ by around 600 time-steps. 
This proved necessary because the third order time integration 
scheme of the Pencil Code is not constructed to conserve orbital 
angular momentum and energ y. Using the highly o ptimiz ed or- 
bital dynamics code SWIFT, ILithwick & ChiangI (l2007h solve 
the same problem with slightly less than five time-steps per or- 
bit. 

In Fig. [21 we show the eccentricity evolution of the particle 
ring. For a coefficient of restitution of 6 = 0.3 the particles relax 
to an equilibrium eccentricity of around ^rms = 0.001, compara- 
ble to dxir. A higher coefficient of restit ution of 6 = 0.6 leads in- 
stead to catastrophic heating of the disc (iGoldreich & Tremaind 
Il978h . with an eccentricity that evolves linearly with time. The 
results presented in Fig. [2 show that the superparticle collision 
algorithm is in excellent agreement with iLithwick & Chiang 
(l2007l) in the limit of inflated particles. 



3.1. Density evolution 

The width of a particle ring increases due to collisional vis- 
cosity. Since the collision time-scale scales inversely with par- 
ticle density, the collisional evolution slows down with time. 
An analytical s olution to the diffusion pro blem was found by 
iPetit & Henonl (Il987h . In the notation of ILithwick & ChiangI 
(12007^ the width cr^ of an initially narrow ring increases accord- 
ing to 



/ 36 {6xf t \ 



1/3 



(11) 



Here ky is a dimensionless factor that depends on the coefficient 
of restitution 6, 6x is the grid spacing, r is the mean radial coor- 
dinate of the p articles, A^tp is the particle n umber and t the time. 

We follow Oth wick & ChiangI (l2007l) and define an initially 
very narrow ring of radial extent 2 A = 10"^. The units follow 
from our choice of GM = 1 . The evolution of the radial width 
is shown in Fig. [3] over 10"^ orbits. We ov erplot the analytical 
solutio n for = 0.016, similar to the fit in ILithwick & ChiangI 
(l2007l) . and find excellent agreement. 




0.000 



2000 4000 6000 
t/yr 



8000 10000 



Fig. 2. The eccentricity evolution of particles orbiting a central 
gravity with GM = 1. Relatively inelastic collisions, with co- 
efficient of restitution e = 0.3, evolve towards an equilibrium 
eccentricity of 10"^, with orbital excursions comparable to the 
grid spacing. More elastic collisions, with 6 = 0.6, lead to catas- 
trophic heating of the particle system . The results follow closely 
Fig. 1 of ILithwick & ChiangI (l2007l) . 



0.010 




0.001 



Fig. 3. The width of a particle ring orbiting a central gravitating 
mass versus time. The 10000 particles were initially placed in a 
ring centred at r = 1 and a width of 2 A = 10 "^, similar to the grid 
spacing. Compare to upper panel of Fig. 3 in ILithwick & ChiangI 
(i2007i) . 



3.2. Superparticle collisions in the local frame 

Hill's equations describe motion relative to a frame that coro tales 
with the Keplerian frequency Q at an arbitrary distance from the 
central gravity source. The coordinate axes are defined such that 
X points radially outwards and y points along the flow of the disc. 
The 2-D equations of motion of particles are 



-\-2Qvy -h 3Q^x, 



dt 

dVy 

— = —2Qvr . 
dt 



(12) 
(13) 



Particle positions are evolved through x = v. The boundary con- 
ditions are periodic in the azimuthal direction. Particles passing 
over the inner (outer) radial boundary get the velocity (3/2)QLx 
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Energy accounted for: 




(0.332736=100%) 

Jacobi constant 

Energy released by stresses 



Energy dissipated by 

inelastic collisions 



10' 



Fig. 4. Evolution of energy in a shearing box simulation where 
particles have a mean-free-path of A = 0.1// and coefficient of 
restitution 6 = 0.3. Drag forces are ignored. The Jacobi constant 
falls due to dissipative collisions. By monitoring the energy re- 
leased as particles pass the boundaries and the energy dissipa- 
tion by inelastic collisions we can account for all the energy in 
the system. 



subtracted (added) to their azimuthal velocity. We also refer 
to the frame as the shearing box. We consider a box size of 
Lx = Ly = 0.2 covered by 32^ grid cells and 102,400 particles. 
The conserved energy (Jacobi constant) is 



1 .2 1 .2 
-mx -h -my ■ 

2 1 ' 



-mQ X 



2^2 



(14) 



Elastic collisions re-orient the particles without changing en- 
ergy, and thus convert circular orbits into eccentric ones while 
conserving energy. Ignoring gas, which damps the velocity rela- 
tive to the gas and hence the eccentricity, elastic collisions con- 
serve the Jacobi energy. Fig. [4] shows the energy of particles 
versus time in local frame simulation with inelastic collisions. 
Particles are initialised with random position and velocity vec- 
tors (^v =1). The mean-free-path is /I = 0.1//, giving an initial 
collision time-scale of Tc ~ 0.1. The coefficient of restitution 
is 6 = 0.3. The Jacobi constant falls with time due to the en- 
ergy dissipated by inelastic collisions. At the same time parti- 
cles passing over the radial boundaries release energy from the 
Keplerian shear through their mean Reynolds stress (the code 
tracks and outputs that energy release for each particle passing 
the radial boundary). All energy in the system is accounted for 
in these three reservoirs. 

3.2.1. Shear during collision 

Particle collisions in the shearing box release energy from the 
Keplerian shear into random motion, leading in the absence of 
drag forces either to catastrophic heating (Vrms ^ oo) or to 
an equilibrium with energy dissipation in inelastic collisions 
(Vrms^^-^ where R is the particle radius). Discounting the for- 
mer option, the result of the latter can be artificially exaggerated 
by the numerical scheme because we identify the collision be- 
tween two superparticle swarms with the collision between two 
members of the swarms located at the respective swarm centres. 



10" 
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10« 

10-^ 
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Fig. 5. Evolution of particle rms speed in the shearing box for 
a simulation with normal collisions (KS, blue/black line) and a 
simulation in which the relative Keplerian shear is subtracted 
when determining the collision time- scale and outcome (NS, 
red/gray line). The top panel shows the decay of initially ran- 
dom particle motion due to inelastic collisions {e = 0.3). The 
rms speed can not fall below Vrms ~ (Sx)Q for KS collisions, due 
to the energy release from the Keplerian shear. In the simulation 
with NS collisions, on the other hand, the rms speed continues 
to decay towards zero. In the bottom panel we consider elastic 
collisions (e = 1.0) with zero random motion initially. Energy is 
released from the Keplerian shear. The blue line shows results of 
simulations with NS collisions, rerun from snapshots of the KS 
simulation at various times. The two solutions match increas- 
ingly well when the particle rms speed increases above (Sx)Q. 



In reality collisions would occur between neighbouring particles 
separated by less than their physical diameter. The naive numer- 
ical algorithm would make the system settle for an equilibrium 
where Vrms~(^-^)'^, where Sx is the grid spacing and also the 
typical distance between superparticle centres. This rms speed 
greatly exceeds the desired Vrms ~ R^- In other words, the naive 
collision algorithm will input artificial heating. 
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Table 1. Simulation parameters. 



Run 




NxXNyX N, 


N 




Collisions 


6 






SI64_nocoll 


0.2 X 0.2 X 0.2 


64 X 64 X 64 


300,000 


0.3 






100 




SI64_eL0 


0.2 X 0.2 X 0.2 


64 X 64 X 64 


300,000 


0.3 


KS 


1.0 


100 




SI64_e0.3 


0.2 X 0.2 X 0.2 


64 X 64 X 64 


300,000 


0.3 


KS 


0.3 


100 




SI64_e0.3_NS 


0.2 X 0.2 X 0.2 


64 X 64 X 64 


300,000 


0.3 


NS 


0.3 


100 


52 


SI128_nocoll 


0.2 X 0.2 X 0.2 


128 X 128 X 128 


2,400,000 


0.3 






50 




SI128_eL0 


0.2 X 0.2 X 0.2 


128 X 128 X 128 


2,400,000 


0.3 


KS 


1.0 


50 




SI128_e0.3 


0.2 X 0.2 X 0.2 


128 X 128 X 128 


2,400,000 


0.3 


KS 


0.3 


50 




SI128_e0.3_NS 


0.2 X 0.2 X 0.2 


128 X 128 X 128 


2,400,000 


0.3 


NS 


0.3 


50 


19 



Col. (1): Name of simulation. Col. (2): Box size in scale heights. Col. (3): Resolution. Col. (4): Number of particles. Col. (5): Friction time. Col. 
(6): Collision type. Col. (7): Coefficient of restitution. Col. (8): Simulation time in orbits. Col. (9): Time of starting self-gravity. 



Collisions between particles of radius R ^ Sx can be mod- 
elled by subtracting the Keplerian shear part from the relative 
speed both for determining the collision time- scale and for deter- 
mining the outcome of the collision. Decomposing the azimuthal 
velocity field sls y = Vy -\- vf\ where v^^^ = -(3/2)Qx is the 
Keplerian shear velocity and V3; is the peculiar velocity, we can 
calculate both the collision time- scale a nd ou tcome in terms of 
Vy (together with Vx and v^). iLyra et all (l2009l) applied a similar 
trick to subtract off' the entire (Keplerian plus peculiar) gas veloc- 
ity from the particle velocity. However, two particles moving at 
the same velocity as the local gas do not necessarily avoid colli- 
sions, even if the gas is incompressible, since the particle motion 
is not completely coupled to the gas. Therefore we choose in this 
paper to subtract off' only the Keplerian orbital speed from the 
particle velocity. The dynamical equations of the Pencil Code are 
already formulated relative to the Keplerian shear, so subtracting 
off' the shear is natural to the governing system of equations. 

Collisions relative to the Keplerian shear conserve both the 
total momentum and the momentum relative to the Keplerian 
shear, but the energy in elastic collisions is only conserved rela- 
tive to the Keplerian shear. To see this, consider the kinetic en- 
ergy of two particles, 

E=lm {vl, + [v,i + vf^f] + {vl, + [v,2 + ^^f] . (15) 

Here m is the mass of a superparticle, assumed to be the same 
for both colliders. An elastic collision solved in terms of (Vxi , v^;! , 
Vx2, Vy2) conserves both the sum of the squares of those velocity 
components, as well as the squares of v^^^ and v^^^ (the latter is 
true since the position x is not changed by the collision). The 
diff'erence in energy before and after the collision is therefore 

AE = £'after " ^before = m[Av3;iV^^^^ -h AV3;2V^2^] . (16) 

This result holds also in 3-D. The energy diff'erence is generally 
not zero, even though Av^;! = -Avy2 by momentum conservation, 
since the off'set v^^^ is not the same for the two particles. The 
non-conservation is nevertheless small: the azimuthal velocity 
change in the collision is uncorrected with the Keplerian shear 
velocity, so (Av3;V^^^)box ~ 0- The particle integrator's slight non- 
conservation of Keplerian orbits is not a serious limitation in 
simulations where the dynamics is driven by hydrodynamical 
instabilities and drag forces. The correct relative Keplerian shear 
based on the physical size of the particles can in principle be 
added artificially, to obtain the correct energy release from the 
shear, but this is negligible for 1-10 cm particles considered in 
this paper. 

The total angular momentum of two colliding particles, 
L = mri X vi mr2 x V2 , (17) 



is conserved in the collisions, both with and without Keplerian 
shear in the collision, as long as the force during the collision 
acts along the line connecting the two particles. This is the case 
both with and without Keplerian shear. For equal-mass particles 
we can write the change in the velocity as Avi = -Av2 = c{r2 - 
ri), giving 

AL = mri x Avi -h mr2 x Av2 = . (18) 

The above arguments for energy and angular momentum con- 
servation are generalisible to distinct particle masses as well. 
However, while the Monte Carlo collision scheme in itself is 
fully consistent with distinct particle masses, correct energy 
equipartition among particle sizes can not be obtained with 
equal-mass superparticles (see discussion in Appendix I A. lb . 

In the following we use the abbreviations KS for collisions 
that include Keplerian shear and NS for collisions where the 
Keplerian shear is subtracted off' when determining the colli- 
sion time-scale and outcome. Fig. [5] shows the evolution of the 
particle rms speed in a shearing box simulation. The top panel 
shows the decay of initially random particle motion by inelas- 
tic (6 = 0.3) collisions for KS collisions and for NS collisions. 
KS collisions decay towards Vrms ~ {6x)Q, the random motion 
released by the Keplerian shear in a single collision. NS colli- 
sions on the other hand continue to decay towards zero. In the 
bottom panel of Fig. [5] we start with zero random motion and 
observe how elastic {e = 1.0) KS collisions heat up the system. 
Rerunning the simulation with elastic NS collisions from various 
starting times of the KS simulation shows clearly that the evolu- 
tion of the system is very similar as long as the particle rms speed 
is larger than (Sx)Q. In actual simulations with gas and hydrody- 
namical instabilities driving particle dynamics with characteris- 
tic motion much faster than v ~ (Sx)Q, one can subtract off' the 
Keplerian shear term when determining the time- scale and out- 
come of collisions and still model the correct system, without 
any spurious energy released by bloated particles. 

4. Particle collisions and the streaming instability 

Armed with a collision algorithm for superparticles, we are now 
ready to explore the eff'ect of particle collisions on particle con- 
centration by streaming instabilities and planetesimal formation 
by self-gravity. The streaming instability feeds off' the relative 
(streaming) motion of gas and particles in protoplanetary discs 
and has a character istic length scale comparab le to the sub- 
Keplerian length rjr (lYoudin & Good man, 2005). Here r/ is the 
radial pressure gradient parameter of|Nakagawa et al. (1986) and 
r is the distance to the central star. IJohansen et al.l (l2009h and 
iBai & Ston3 (1201 Obi) demonstrated that the streaming instability 
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Fig. 6. Maximum particle density, relative to the mid-plane gas 
density, versus time for a series of 64^ simulations (top plot) 
and 128^ simulations (bottom plot) of turbulence driven by the 
streaming instability with different treatment of collisions. The 
maximum particle density increases by a factor approximately 2 
when doubling the resolution, but the maximum density peaks 
are consistently 50% lower when including particle collisions. 
Note the different scale of the axes in the two plots. 
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Fig. 7. Maximum particle density, relative to the mid-plane gas 
density, versus time for simulations with normal collisions (KS) 
compared to simulations where we subtract off the Keplerian 
shear difference between particle pairs when calculating the col- 
lision time and the outcome of the collision (NS). NS collisions 
display more than three times higher particle densities than KS 
collisions. Peak concentrations fill a larger fraction of the simu- 
lation time at 128^. 



leads to strong particle clumping when the heavy element abun- 
dance of the disc is above a threshold value of Z ^ 0.02 for par- 
ticle si zes Qzf > 0.1 (and moderate radial drift, see lBai & Stond 
l2010d) . Clumping proceeds as initially very low amplitude parti- 
cle overdensities accelerate the gas towards the Keplerian speed, 
hence reducing the local head- wind, which in turn slows the ra- 
dial drift of the particles. Drifting particles pile up where the 
head- wind is slower, causing exponential growth of the particle 
density as the pa rticles continue t o incr ease their drag force influ- 
ence on the gas. I Johansen et al.l (l2009l) found that overdense re- 
gions contract when including particle self-gravity and that even- 
tually a number of gravitationally bound clumps form. These 
models nevertheless did not include any particle collisions. 

We perform 3-D simulations where the gas is modelled on 
a fixed grid and solid particles with superparticles. We solve the 
standard shearing box equa tions for gas and particles (same as in 
iJohansen & Youdinl l2007l but with additional vertical gravity). 
The frame rotates at the Keplerian frequency at a fixed orbital 
distance r from the star. The coordinate axes are oriented such 
that X points radially outwards, y points along the rotation direc- 
tion of the disc, while z points perpendicular to the disc along £1. 
The gas is subjected to a radial pressure gradient which reduces 



its orbital speed by the positive amount Av = 0.05 Cs. Particles do 
not feel this radial pressure gradient, and the resulting relative 
mo tion between particle s and gas drives the streaming ins tabil- 
ity (iGoodman & Pindoii l200Ql: lYoudin & GoodmanL l2005b . We 
consider a cubic box with side lengths Lx = Ly = = 0.2H, 
where H = c^jQ is the gas scale height, to capture the fastest 
growing modes of the streaming instability of marginally cou- 
pled particles, Asi/H ~ T]r/H ~ Av/Cs = 0.05. This is also the 
characteristic scale of Kelvin-Helmholtz instabil ities, thriving in 
the vertical shear in the gas and particle velocity (Youdin & ShuL 
2002; Lee et al, 2010), although Bai & Stone (2010b) demon- 
strated that the streaming instability is dominant over Kelvin- 
Helmholtz instabilities in setting the dynamics of particle layers 
with-QTf > 0.1. 

The friction time of the particles is fixed at Qzf = 0.3 in 
all simulations, corresponding to approximately 20-cm rocks 
around the locatio n of the asteroid belt a t 3 AU, and to 6-mm 
pebbles at 30 AU (IWeidenschillingL Il977l) . The particle column 
density is set to 2% of the total gas column density, the latter 
including the gas beyond the vertical boundaries of the box. For 
our choice of Av strong particle clumping can only be obtained 
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Fig. 8. The three components of the particle velocity as a func- 
tion of the radial position within a grid cell. Three grid cells were 
chosen at ^ = 45rorb of the run SI128_e0.3, one with the highest 
particle density in the box, one with a particle density close to 
100 times the gas density and finally one with a particle density 
close to 10 times the gas density. Both systematic and random 
particle motion is present within the grid cells. The Keplerian 
shear is clearly visible in the velocity (marked with a solid line 
in the middle panel). The cells with the highest density have gen- 
erally a slower random motion and are thus more aff'ected by the 
Keplerian shear. 



at such super- solar metallicit}0. The average dust-to-gas ratio in 
a box of L, = 0.2// is <Pp/pg) ~ 0.25 when Z = 0.02. We set 
sound speed Cs, Keplerian frequency Q and mid-plane gas den- 
sity po to unity, so these form the natural units of the simulations. 

We compare results obtained without and with particle colli- 
sions. Simulations with particle collisions are run in three varia- 
tions: either with elastic collisions {e - 1.0), with inelastic colli- 
sions (e - 0.3) or with inelastic collisions where Keplerian shear 
is subtracted off' when determining the time- scale and outcome 
of collisions. Simulation parameters are given Table [T] Each par- 
ticle swarm contains a mass per volume of Pp/po ~ 0.219 for the 
considered particle number at both 64^ and 128^. 



^ The threshold for clumpin g can be estimated analyt ically to be Z - 
rfrIK ) (lYoudin & S hu. 2002!). lBai &"St one (20103) and lJohansen et all 
(l2007h confirmed numerically that the threshold for particle clumping 
by the streaming instability shifts towards higher (lower) metallicity as 
the sub-Keplerian speed difference Av is increased (decreased). 



4A. Maximum particle density 

We monitor the maximum particle density regularly in the sim- 
ulations. In Fig. [6] we show the maximum particle density ver- 
sus time in simulations with 64^ grid cells and 128^ grid cells, 
respectively. Simulations without collisions generally achieve 
higher particle density - up to 600 times the gas density at 64^ 
and 1200 times the gas density at 128^. Elastic collisions and 
inelastic collisions with 6 = 0.3 give very high particle densi- 
ties too, but the peaks have an approximately 50% lower value 
than in simulations without collisions. Elastic collisions achieve 
a somewhat lower maximum density than inelastic collisions. 
The kinetic energy dissipation in inelastic collisions reduces the 
random motion of the particles and allows higher particle con- 
traction. 

The inclusion of Keplerian shear during the collision can 
lead to unphysical results, since the shear term is exaggerated 
by enlarging particles to the size of a grid cell. The exaggerated 
kinetic energy input will in turn suppress concentration peaks, 
in agreement with what is seen in Fig. [6l In Fig. [7] we show the 
maximum density in simulations with inelastic KS and NS col- 
lisions respectively (and the results without collisions for com- 
parison). Simulations with NS collisions display a three times 
higher maximum density than simulations with KS collisions. 
The maximum density is even a factor 2-3 times higher than in 
simulations without collisions. This way collisions actually pro- 
mote particle concentration. 

In Fig. [8] we analyse the particle motion within three grid 
cells of the run SI_128_e0.3. We choose the grid cell with the 
maximum particle density in the box and two grid cells with a 
particle density close to 100 and 10 times the gas density, respec- 
tively. The particle velocity shows both systematic trends and 
random motion within the cells. The random motion is slower in 
the cells of higher density. The Keplerian shear is clearly visible 
in the velocity of particles in the two densest grid cells. Thus 
the hydrodynamical simulations are prone to spurious heating, as 
explained above. Subtracting off' the Keplerian shear term when 
determining the time- scale and outcome of collisions avoids this 
problem. Fig. [8] also shows a systematic trend in the radial par- 
ticle velocity. Radial convergence and divergence in the particle 
velocity are expected when particles concentrate in radial bands 
and when the concentrations dissolve again. We do not attempt 
to correct for this systematic velocity within grid cells, but note 
that systematic trends from smooth gradients will decrease with 
increasing resolution. 

4.2. Particle concentration versus scale 

Overdense particle sheets co ntract radially under the action of 
self-grav ity and dr ag forces (lYoudini 120111: iMichikoshi et al.L 
[iOlO; ShariflL&Cuzzi, 20li). A full non-axisymmetric collapse 
is initiated when the particle density crosses the Roche density 
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(19) 



The mass of the planetesimal will be characterized by the scale 
over which the Roche density is achieved. To quantify the scale- 
dependence of the particle concentrations, we measure the max- 
imum particle density over cubic regions of side length A^t grid 
cells, increasing A^t from 1 to A^;c- We ensure that all concen- 
trations centres are probed by stepping the measurement re- 
gion through the entire grid. Measurement regions crossing the 
boundaries are handled by expanding the particle density field 
with its periodic counterpart in all directions (glueing together 
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Fig. 9. Maximum particle density, relative to the mid-plane gas 
density, as a function of scale, for simulations with NS collisions 
(top panel) and simulations with no collisions (bottom panel). 
Diamonds indicate the maximum density over a given scale, 
while pluses indicate the mean of the time-dependent maximum 
density. Simulations with NS collisions display good conver- 
gence in the maximum density, following closely a max(pp) oc 
law (thin black line), while the mean of the maximum den- 
sity increases from 64^ to 128^, due to a higher temporal fill- 
ing factor of major concentration events at higher resolution 
(see Fig. O. The dashed line shows the maximum density for 
a uniform razor- thin mid-plane layer for comparison. Blue dot- 
ted lines show the Roche density for the minimum mass solar 
nebula at 3 AU from the central star, and for five times less and 
more massive nebulae. The red dotted line indicates the charac- 
teristic length scale of the streaming instability, L = T]r. Particle 
densities above 10^ times the gas density are reached in regions 
smaller than ^ 0.003/f, equivalent of L ^ 50,000 km at 3 AU. 



3^ copies which are identical except for a shift due to Keplerian 
shear). 
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Fig. 10. Zoom in on the densest grid cell in SI128_e0.3 JMS at 
t = 32rorb. The overdense particle structure is elongated along 
the shear direction with a density decreasing in all directions 
from the densest point. The lower-right panel shows the particle 
density average over shells of thickness one grid cell and a 1/r^ 
power-law overplotted. 



For snapshots saved once per orbit from t = 20rorb to 
t = SOJorb we calculate the maximum particle density as a 
function of scale. The results are shown in Fig. [9] for simula- 
tions with NS colHsions (SI64_e0.3JS[S and SI128_e0.3JS[S) in 
the top panel and simulations with no collisions (SI64_nocoll 
and SI128_nocoll) in the bottom panel. We extend the measure- 
ments of SI64_e0.3 JMS to ^ = 60rorb to catch a major concen- 
tration event (see top panel of Fig.l?]). We indicate in Fig.[9]both 
the maximum density over all times and the mean of the time- 
dependent maximum density. The maximum scale-dependent 
density in NS simulations is very similar at 64^ and at 128^. This 
quantity is nevertheless very sensitive to the low-number statis- 
tics of the concentration events. A more robust measure is the 
mean of the maximum density. This measure increases some- 
what from 64^ to 128^. It is also evident from Fig. [7] that ma- 
jor concentration events have a higher temporal filling factor at 
128^. Whether this is intrinsic to the streaming instability dy- 
namics or just an eff'ect of running simulations for too short time 
is not possible to discern. 

The apparent linear decrease of logarithmic density with log- 
arithmic scale implies max(pp) oc L"^ as a good model for 
the scale-dependence of the maximum density. Two limits can 
immediately be put on a. The lowest value would stem from 
a razor-thin particle mid-plane layer of uniform density, with 
M oc l}^ giving max(pp) oc MjL^ oc and thus a = \. 
Concentration of all particles in a single point would yield the 
upper limit of a = 3. We overplot in Fig. [9] with a thin black line 
the power law max(pp) oc L"^, fitted to match the mean density 
of the box at L = 0.2//. The a = 2 power law follows the data 
extremely well. This implies that M oc L, i.e. that the particles 
primarily concentrate either in 1-D filaments or in spherically 
symmetric clouds of density p(r) oc known i n star forma- 
tion as the singular isothermal sphere solution (e.g. lShulll977l) . 
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In Fig.[TOlwe show the particle density around the densest grid 
point in SI128_e0.3JS[S at ^ = 32rorb. The overdense structure 
appears elongated along the j-direction with the density falling 
rapidly towards all directions (although slower along y). 

Simulations without collisions (bottom panel of Fig. O show 
similar trends as the simulations with NS collisions, but there is 
a marked decrease in the maximum density over the smallest 
shared scale between 64^ and 128^. Nevertheless the mean of 
the maximum density agrees between the two resolutions. 

The convergence in scale-dependent maximum density 
shows that the dynamics of the streaming instability concen- 
tration events is well-resolved and independent of dissipation 
scale and viscosity. This is in contrast to turbulent concentra- 
tion in driven isotropic turbulence which, for a given parti- 
cle size, appears on length scale s that are f ixed relative to the 
Kolmogorov (viscous) scale (Ho gan & Cuz zi, 2007; Pan et al.L 
I2OIII) . In contrast the streaming instability is fixed relative to 
the sub-Keplerian scale r/r ~ 0.05H. At ^ ~ 0.00 16if, probed 
only at 128^, the maximum density in simulations with NS col- 
lisions reaches more than three thousand times the gas density. 
Higher resolution simulations will be needed to test if the par- 
ticle density continues to follow the max(pp) oc L"^ trend, or 
eventually fin ds a smallest s cale. T he 2-D streaming instability 
simulations o f lBai&Stond (l2010ah converged in density statis- 
tics at between 512^ and 1024^ grid cells. Reaching those reso- 
lutions in 3-D is very computationally demanding, but should be 
an important priority for the future. 

5. Planetesimal formation 

The gravitational potential field of the particles is found by map- 
ping the particle density on the grid, using a second order spline 
interpolation scheme, and sol ving the Poisson equa tion using a 
fast Fourier transform method (iJohansen et '^ 120071) . The gravi- 
tational acceleration is interpolated back to the particle positions 
using second order spline interpolation. The strength of the grav- 
ity is defined by the non-dimensional parameter 

which is relat ed to the thin-disc self-gravity pa rameter Q through 
O ^ 1.6G-^ (ISafronovi [19601: iToomreL 11964 . The solar nebula 
of lHavashil (Il98lh has G ^ 0.04 at 3 AU from the sun, the pa- 
rameter depending weakly on the distance. We use G = 0.1 as a 
reference choice in the simulations, but experiment with G down 
to 0.02. 

The total particle mass in the box is 

Mp = <pp)L^ ^ O.OO2//V0 , (21) 

where the mass unit Mq = H^po depends on the temperature and 
location in the disc [H] and the strength of the self-gravity [po = 
(47iG)~^GQ^]. While the expression in Eq. (|2T1) does not depend 
on G, in units where // = po = 1, the p hysical niass un it does. In 
a nebula with the scale-height given bv lHavashil(ll981h . we have 
at r = 3 AU with G = 0.1 a mass unit of Mq ~ 1.3 x 10^^ g and 

Mp ^ 2.8Mceres. 

We activate particle self-gravity in simulations of the stream- 
ing instability with inelastic NS collisions, at times when there 
is little particle concentration, to catch the simultaneous action 
of streaming instability and self-gravity during the next con- 
centration event. In SI64_e0.3_NS we thus start self-gravity at 
t = 52rorb, while in SI128_e0.3JMS we start self-gravity at 



t = 19rorb (see Fig. [7]). We then evolve the simulation for an- 
other 5 orbits, either ignoring collisions or applying the usual 
variation of collision types (elastic, inelastic KS, inelastic NS). 

Results of 64^ simulations are shown in Fig. [TT] Between 
3 and 4 clump^ initially condense out of the dominantly ax- 
isymmetric filament forming by the streaming instability. These 
clumps have masses between a tenth and a third of the dwarf 
planet Ceres - corresponding to contracted radii between 220 
and 330 km, assuming an internal density of 2 g/cm^. All the 
clumps form in a single planetesimal-formation event shortly af- 
ter the onset of self-gravity. The clumps continue to grow mainly 
by accreting particles from the turbulent flow, but no new grav- 
itationally bound clumps form. Clumps eventually collide and 
merge in all simulations. Such clump merging is likely an un- 
physical efl'ect driven by the large sizes of the planetesimals. The 
self-gravity solver does not allow gravitational structures to be- 
come smaller than a grid cell, and that leads to artificially large 
collisional cross sections. A more probable outcome of the real 
physical system is gravitationa l scattering and/or formation of 
binaries (iNesvorny et al .1120101) . 

Results at 128^ are shown in Fig. [121 At higher resolution the 
number of clumps condensing out is about twice as high com- 
pared to the lower resolution simulation. However, the masses 
of the most massive clumps are very similar to lower resolution 
(although a bit higher - up to 60% of Ceres), so it appears that 
higher resolution simply allows lower-mass clumps to condense 
out as well. The masses of the clumps condensing out at 128^ 
resolution correspond to contracted radii between 84 and 405 
km. The ability to form smaller clumps at higher resolution is ex- 
pected from the picture that a radial contraction phase is needed 
before the Roche density can be achievecQ. Higher resolution al- 
lows contraction to narrower bands and thus formation of less 
massive planetesimals. It is nevertheless difficult to compare the 
planetesimal masses condensing out at the two resolutions as the 
initial conditions are not the same. 

iRein et aP (l2010h observed in their 2-D shearing sheet simu- 
lations that inclusion of collisions would lead to condensation of 
fewer and more massive clumps, when compared to simulations 
without collisions. Our Fig. [12] also shows that the simulation 
with no collisions makes the highest number of clumps of all 
the four simulations. Nevertheless the characteristic mass of the 
most massive clumps appears indifferent to the treatment of col- 
lisions. 

Since G controls the relative strength of self-gravity, results 
obtained with a given G can not be scaled to other values of 
G. We vary the self-gravity parameter in 128^ simulations in 
Fig. \T3\ starting self-gravity at the same time as in Fig. [T2l 
Weaker self-gravity gives lower clump masses, but gravitation- 
ally bound clumps of up to 0.01 Ceres masses (or 100 k m radius) 
conde nse even at G = 0.02. The solar nebula model of lHavashil 
(119811) has G ^ 0.04 at 3 AU from the sun. Thus the streaming 
instability allows planetesimal formation in disc models that are 
similar in mass to the solar nebula, in contrast to recent simu- 
lations of planetesimal formation in pressure bumps excited by 
the magnetorotational instab ility which required d isc masses up 
to 10 times the solar nebula (I Johansen et al.Ll20Tl1) . 



^ The algorithm for identifying bound clumps is based o n 2-D col- 
umn d ensity snapshots and is described in detail in Joha nsen et aLl 

dmll). 

^ A similar order of events is seen in simulations of star formation 
in self-gravit ating accretion d i scs aro und supermassive black holes, see 
e.g. Fig. 3 of I Alexander et all (120081) . 
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Fig. 11. Particle column density versus time after self-gravity is turned on at to = 52rorb = in the simulation 

SI64_e0.3 JS[S. An overdense sheet forms by the streaming instability and breaks up in a number of gravitationally bound clumps. 
We indicate the number of clumps and their masses, in units of the mass of the dwarf planet Ceres, in the lower left part of the plots. 
Between 3 and 4 clumps condense out independently of how collisions are treated, with masses slightly smaller than Ceres. Clump 
merging, likely driven by the artificially large sizes of the planetesimals, reduces the number of clumps with time in all cases. Note 
that the initial condition for all four simulations is taken from SI64_e0.3_NS. 



The presented simulations do not catch the transition from 
bound clump to solid planetesimal. However, 'Ne svorny et al.l 
(|201Q|) simulated the gravitational collapse of spherical particle 
clouds and generally found formation of binary planetesimals, 
with the two largest bodies containing a significant fraction of 
the mass of the cloud. The fact that the masses of the most mas- 



sive bound clumps in our simulations are relatively independent 
of resolution allows us to critically compare the mass distribu- 
tion of the clumps to to the observed properties of the asteroid 
and Kuiper belts and extrasolar debris discs. 
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Fig. 12. Same as Fig.[TTl but for 128^ simulations with self-gravity started at = 19rorb = 1 19.381-Q"^ More clumps form initially, 
but the most massive clumps have similar masses to the 64^ simulation. The run with no collisions forms more low-mass clumps 
than the other runs. The initial condition for all four simulations is taken from SI128_e0.3_NS. The total particle mass in the box is 
approximately 2.8 Ceres masses. 



5.1. Application to the Kuiper belt 



The physical mass of the clumps depends on location in the disc 
and on the self-gravity parameter G. While the simulations are 
dimensionless, the translation to physical mass involves multi- 
plication by the mass unit Mq = p^H^ = GQ^H^ /(47iG). In a 
nebula with constant G and T oc r"i/^, the mass unit scales as 



Mo oc r^/"^, so re-scaling to the Kuiper belfl gives planetesimal 
masses 5-6 times higher than in Fig.[TT]and Fig. [121 Contracted 
radii at the location of the Kuiper belt are approximately 80% 
higher than in the asteroid belt, yielding planetesimal radii be- 

^ The orbits of trans-Neptunian objects extend to several 10 AU be- 
yond the orbit of Neptune, although many of these must have formed 
within the orbit of Neptune and been scattered outwards later. Thus we 
take 30 AU as an approximate distance scale. 
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Fig. 13. Evolution of maximum planetesimal mass (full line) and 
total mass in planetesimals (dash-dotted line) for 128^ simula- 
tions with inelastic NS collisions (thin yellow line shows the 
G = 0.1 simulation without collisions for comparison). Colors 
indicate G = 0.02,0.05,0.1. Extended wiggles in the G = 0.1 
curve arise during clump merging. The total particle mass in the 
box is 2.8, 1.4 and 0.56 Ceres masses, in order of decreasing G. 



tween 150 and 730 km. The upper range is com parable to the 
masses of the large st known Kuiper belt objects (IChiang et al.L 
120071: lBrownLl2008h . 

This extrapolation is only valid for an assumed constant 
self-gravity parameter G. The minimum mass solar nebula, with 
Z oc r"^/^, has G oc r^^"^. The weak dependence on radial distance 
from the star gives in the Kuiper belt at r = 30 AU a 10^ 1.8 
times larger G than in the asteroid belt. From Fig. [13] we read off 
an approximate doubling in planetesimal mass when increasing 
G from 0.05 to 0.1. We expect that this scaling holds for larger 
G as well. This way the minimum mass solar nebula gives some- 
what higher masses in the Kuiper belt compared to the constant- 
G extrapolation presented above. 

The comparison to observed planetesimal belts is never- 
theless complicated by a potentially very efficient accretion 
of unbound particles (pebbles and rocks) by the newly born 
planetesimals after th eir formation (I Johansen & Lacerda 120101: 
lOrmel & Klahiil20TQl) . an epoch not captured in our simulations. 
It is interesting to note that, given the power of the streaming 
instability in producing Ceres-mass planetesimals from pebbles 
and rocks, the challenging question may not be how these plan- 
etesimal belts fornfl or how the characteristic mass arises, but 
rather why the planetesimals did not immediately continue to 
grow towards terrestrial planets, super-Earths, and cores of ice 
and gas giants. Perhaps these planetesimal bursts were "aban- 
doned" by the particle overdensity from which they formed, by 
radial drift of the particles, stranding as planetesimal belts. Such 
stranding is evident in the last frames of Fig. [TT] and Fig. [12] 
where the gravitationally bound clumps clearly lag behind the 



^ This does require sufficient amounts of pebbles and rocks to 
begin with, the form ation of which is not yet well-understood 
(lBlum&Wurml . [2008h . 



overdense particle filament. The lag might have be even more 
pronounced if the particle clumps would not be bloated to fill a 
grid cell. 

This stranding scenario is an alternative to the more classical 
view that the asteroid and K uiper belts were disturbed by the 
presence of giant planets (e.g. lKenyon & Bromleyl[2QQ4l) . 

6. Summary and discussion 

This paper focuses on the effect of momentum exchange and 
energy dissipation in collisions on particle concentration by the 
streaming instability and on the subsequent gravitational col- 
lapse to form dense clumps and planetesimals. We develop a new 
algorithm for tracking collisions between superparticles repre- 
senting swarms of physical particles. The time-scale for a par- 
ticle in a given swarm to collide with a particle from another 
swarm is calculated for all superparticle pairs in a grid cell. 
Collisions occur instantaneously if a random number is less than 
the ratio of the simulation time-step to the collisional time-scale, 
ensuring that superparticles collide statistically on the correct 
time- scale. We have demonstrated that this algorithm can be in- 
corporated into a hydrodynamical code at a modest computa- 
tional cost. This is true even for large particle numbers, since 
the number of possible collision partners that are considered in a 
given timestep can be reduced with little or no loss of generality. 

Collisions can have a number of effects on particle dynam- 
ics, by making particle motion more isotropic and by dissipative 
collisions which drain kinetic energy from the system. We have 
considered the simplest case of a constant coefficient of restitu- 
tion (either unity or 0.3), but a more physically motivated coeffi- 
cient of restitution, depending on material properties and impact 
speed and angle, could be easily implemented in the scheme. We 
emphasize that we have focused in this paper entirely on parti- 
cles with a friction time of 0.3 relative to the local Keplerian 
time-scale, corresponding to 20-cm rocks in the asteroid belt 
and 6-mm pebbles at 30 AU. Future studies will be needed to 
determine the influence of particle collisions on the dynamics of 
smaller and larger particles and on their ability to form planetes- 
imals. 

Our simulations show that collisions are important to con- 
sider when modelling particle concentration by the streaming 
instability. Taking into account the energy dissipation in inelas- 
tic collisions increases the maximum particle density. This in- 
crease is most pronounced, more than a factor of three compared 
to simulations with no collisions, when we ignore the relative 
Keplerian shear for determining the collision time- scales and 
outcomes. We argue that the Keplerian shear velocity should 
be subtracted when determining the outcome of collisions be- 
tween superparticles representing physical particles that are 
much smaller than a grid cell. The collision algorithm enlarges 
particles to the size of a grid cell during a collision, and this 
can lead to unphysical heating of the particle component if the 
Keplerian shear is included during the collision. 

The treatment of collisions has no apparent effect on the 
planetesimals which form by self-gravity. The masses of the 
most massive planetesimals are relatively independent of the in- 
clusion or absence of collisions, although we find some evidence 
that more low-mass clumps condense out in simulations without 
collisions. The particle densities reach several hundred and even 
thousand times the gas density both with and without collisions - 
much higher than the Roche density which governs gravitational 
collapse - and that may explain why particle collisions play a 
relatively small role in determining the outcome of the gravita- 
tional contraction to form planetesimals. The simulations show a 
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characteristic planetesimal mass-scale comparable to the dwarf 
planet Ceres at the location of the asteroid belt. The mass-scale 
increases approximately linearly with distance from the central 
star, giving almost double the contracted radius at the distance of 
the Kuiper belt. This scaling may explain why the largest Kuiper 
belt objects are bigger than the largest asteroids. 

Particle collisions are also important as a stepping 
stone towards implementing coagul ation and fragmentation 
in planetesimal formation models (lOrmel & Spaansl 120081: 
IZsom & Dullemondi l2008l) . Including all the physics relevant 
for modelling particle-dominated self-gravitating flows is a ma- 
jor task, but the reward will be a much better understanding of 
the important step from pebbles and rocks to planetesimals and 
dwarf planets. 
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Appendix A: Collision time from friction time 

In connection with the presence of gas it is convenient to express 
the collision time-scale in terms of the gas friction time-scale. In 
the Epstein drag force regime, valid when the radius of a particle 
R is smaller than ( 9/4 tim es) the mean free path of gas molecules 
(Weidens chilling.. 19771) . the friction time-scale is 



Tf 



Rp. 

CsPg 



(A.l) 



Here p, is the material density of the particles, while Cs and pg 
are the sound speed and density of the gas molecules. 

The time- scale for a particle of radius Rk to collide with a 
swarm of particles with physical radius Rj is 



1 



njCTjkSvjk 



(A.2) 



where o-jk is the mutual collisional cross section. Writing fur- 
ther hj = Pj/rrij and ajk = 7i(Rj -h Rk)^ and assuming spherical 
particles we arrive at 



(4/3)p.7?5 



' Pj{Rj + Rkfdvjk' 
In terms of the friction time we get 



Jk) _ 1^0)^ _^ 



Si) 



(A.3) 



(A.4) 



f / 



For collisions between equal-sized particles, with r^^^ = Tf \ the 
expression reduces to 



3 pj Svjk 



(A.5) 



A time-dependent numerical solution of a collisional particle 
system must take collisions into account when choosing the 
time-step. The time-step criterion of the Monte Carlo collision 
scheme originates in the requirement that two particles can col- 
lide at most once during a time-step, i.e. the collision probability 
P = St/Tc between any two particles in the same grid cell must 
be much smaller the unity. This time- step is independent of the 
maximum density in a grid cell, since particles in dense grid 
cells have many collision partners and hence can suffer more 
collisions in the same time- step. 
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In the streaming instability simulations presented in Sect. IH 
and Sect. Owe observe a typical particle rms speed Sv ~ 0.025cs. 
The mass density represented by a single superparticle is pp ^ 
0.219pg and the friction time is Qrf = 0.3 (we normalise here 
by the Keplerian frequency Q which we define in Sect.©. This 
gives Qtc ~ 18 from Eq. (IA.5I) . The Courant criterion for the 
hydrodynamical part of the streaming instability gives the time- 
step ^%dro = 0.000625^2-1 for 64^ and ^^hydro = 0.0003 125i2-i 
for 128-^ simulations. Therefore we can ignore the collision time- 
scale in the simulations when determining the numerical time- 
step. 



A.1. Multiple particle sizes 



Eq. (IA.2I) defines the collisional time- scale between particles of 
two sizes. For two superparticles of equal internal particle num- 



ber (h) we have rf^ = r^^ because the cross section cTjk and 
relative speed Svjk are symmetric in (j, k). However, equal par- 
ticle number per superparticle is numerically expensive, as the 
mass of a superparticle in that case scales as R^, requiring many 
more superparticles to represent an equal mass of smaller par- 
ticles. The second complication is that the collision time-scale 
becomes very short for smaller particles. 

A more common approach is to have equal mass per super- 
particle. In that case we can define a collision time- scale as the 
time for all mass in particle j to interact with all mass in particle 
k. This time- scale is shared between the two particle species and 
is given by 



-max(T, 



(A.6) 



To illustrate this, take small particles of friction time r^^ 

and large particles of friction time T^^^ = 100. The collision time- 
scale for the large particles is 100 times shorter than for the small 
particles, because the superparticle with small physical particles 
contains 100 times more particles in the swarm. However, the 
time- scale for collision between a large and a small particle does 
not imply that all small particles have collided during that time. 
The correct time- scale is the time for small particles to collide 
with large particles. When an average small particle has expe- 
rienced a collision, then all small particles have collided with 
a large particle, and all the mass in the two superparticles have 
interacted. 

After waiting the common collision time Tc, the collision 
outcome can be solved as if the two colliding particles had equal 
mass, since eff'ectively a large particle collides with mk/mj small 
particles during this time. This approach is slightly inconsistent 
since discrete collisions with N particles of mass mj is not equal 
to a single collision with a particle of mass Nrrij. A collision be- 
tween a particle of velocity Vk and a stationary particle results in 
the new velocity 
mk-mi 

v[ = — -Vk. (A.7) 

mk-\- rrij 

After N such collisions the velocity of particle k is 

(\N 
-] Vk. (A.8) 

In the limit where Vk-v'j^ = Avk «c Vk, this equation describes a 
velocity damping 

— = ■ Vk (A.9) 

Tc ruk + rrij 



with characteristic time-scale Td = Tdnik + mj)/(2mj). 
Completely braking down the large particle requires infinite 
time, whereas a single discrete collision with an equal-mass par- 
ticle would remove all the momentum from particle j in one col- 
lision time. 

To really get the collisional energy equipartition right be- 
tween particles of diff'erent sizes one would have to allow for 
collisions between a large particle and individual smaller parti- 
cles. This could either be done by letting superparticles not rep- 
resent the same mass, but rather the same number of particles. 
However, this approach would become unpractical to model a 
large span in particle sizes, since a huge number of superparti- 
cles would be required to represent the low-mass particle bins. 
Alternatively the collision between a swarm of large and small 
particles could be modelled on the collision time- scale of indi- 
vidual collisions, distributing afterwards the energy and momen- 
tum of the particle that suff'ered the collision among the entire 
swarm of small particles or among all particle swarms within 
its mean free path. However, for a large span in particle sizes, 
this would still require a very small time- step and is therefore 
unpractical. We simply note here that while collisions between 
unequal- sized particles can be modelled with the right conser- 
vation properties, actual equipartition of particle energies would 
require an adaptation of the collision algorithm. 

Appendix B: Limiting the collision number 

During the gravitational contraction of particle clumps the num- 
ber of particles in a grid cell can become very large, on the order 
of 1000s of 10000s. Tracking (1/2)A^(A^ - 1) possible colHsions 
per grid cell then becomes very computationally expensive. 

However, particles do not collide with all possible partners 
during a single time- step. One can limit the number of collision 
partners, while maintaining the overall collision rate, by sam- 
pling only a subset //max of the possible collisions. Considering 
only A^max out of the - 1 collision partners for each particle 
in a grid cell, while increasing the collision probability for each 
collision partner by (A^ - l)/A/max, yields statistically the same 
number of collisions. 

Consider as an example 101 particles in a grid cell, with the 
collision probability between any two particles of 10"^. Particle 
1 will then on the average collide with 1 other particle. However, 
calculating the collision probability with 100 other particles is 
expensive, even when it does not lead to a collision, which is 
most often the case. Instead we let particle 1 only interact with 
particles 2 to 6, and give each collision the probability 10"^ in- 
stead of 10"^. Particle / has particles / 1 to / -h 5 as collision 
partners. When reaching particle 97, the collision partners wrap 
around to particle 1 again, and this way all particles on the av- 
erage get 10 collision partners (5 of higher index and 5 of lower 
index) instead of 100. 

When reducing the number of collision partners, one has to 
be careful that the particles do not interact only with particles 
of a nearby index in each time- step. To avoid any such spurious 
particle preferences, we therefore shuffle the order of particles 
inside a grid cell in each time- step. We have empirically found 
that reducing the number of collision partners becomes impor- 
tant when there are more than 100 particles in a grid cell. We 
show in Fig. IB. II the rms speed of particles undergoing inelas- 
tic collisions with coefficient of restitution e = 0.3. We use 100 
particles per grid cell and show results where we consider all 
particles in a cell to be collision partners together with results 
where we limit the collision partners to 10 and 2. The results 
are indistinguishable, but the code speed is significantly higher 
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Fig.B.l. Evolution of particle rms speed in simulation start- 
ing with random motion of amplitude 1. Particles have mean- 
free-path /I = 0.1 and coefficient of restitution e = 0.3. Drag 
forces are ignored. The blue line shows the results of a simu- 
lation with 100 particles per grid cell and full collision partner 
list, while the red and golden lines show the results of limiting 
the collision partners to 10 and 2, respectively, while increasing 
the collision probability accordingly. The results are extremely 
similar. The lower panel shows the instantaneous inverse code 
speed. Limiting the number of collision partners has increased 
the speed by a factor of approximately three. 



when limiting the number of collision partners (lower panel of 
Fig. IB. lb . The typical speed of the Pencil Code for a hydrody- 
namical simulation with two-way drag forces between gas and 
particles is ~10 yus per particle per time-step. Fig. IB. II shows 
that the computational time needed for superparticle collisions 
is similar to or lower than the time needed for gas hydrodynam- 
ics, particle dynamics, and drag forces, if the number of collision 
partners is kept below approximately 100. 

A side effect of reducing the number of collision partners is 
that the maximum number of collisions is reduced accordingly. 
Therefore it be must required that the boosted collision proba- 
bility P' = P(N - 1)/A^max is always much smaller than unity. 
This must hold for all particle pairs. One can use the maximum 
relative speed between any two particles within a grid cell to 
estimate the smallest allowed A^max that keeps all P' ^ I. 

Each swarm in our simulations represents Pp/pg ^ 0.219. 
The base probability for collision between two superparticle 
swarms with random motion Sv/c^ ~ 0.025 is P = St/Zc ~ 10"^ 
using Eq. (IA.5I) and a typical hydrodynamical time- step St at 64^ 
and 128^ resolution. The maximum density reached in the sim- 
ulations is Pp/Pg ~ 3000 (see Fig.l?]), giving ^ 13700 particles 
in the densest cells. We use A^max = 100 and thus the maximum 
boosted probability is P' ~ 10"^, safely in the regime where the 
collision time- scale can be ignored when determining the nu- 
merical time- step of the code. 



Appendix C: From superparticles to inflated 
particles 

Consider a particle component of mass density pp. A superparti- 
cle can maximally hold a particle number N (equivalently parti- 
cle mass density pp) that covers the whole area of the grid cell. 



AT ^3 Pp (^^y 

No- = — (Sx)o- = 

mp Pp A 



< (6xr 



(C.l) 



Here cr is the cross section of a swarm member and A is the 
mean free path of physical particles in the system. This gives a 
maximum superparticle mass density of 



Pp =Pp 



Sx 



(C.2) 



At this mass density the Monte Carlo method breaks down be- 
cause the superparticle area is larger than a single grid cell (this 
is not taken into account in the model because collisions are only 
considered when superparticles share the same grid cell). The 
free path of a test particle encountering this maximum density 
superparticle is 



A ■ 



1 

ho- 



Pp 
Pp 



A = Sx, 



(C.3) 



using (T = l/(An) = l/(An). Thus the maximum area criterion 
coincides with the particle density where the free path is the 
same as the grid spacing, giving a collision probability of unity 
when the particle enters a grid cell occupied by a superparticle. 
This is in fact equivalent to the inflated particle approach, i.e. 
that overlapping particles always collide. 

We still must show that the mean free path of the system is 
equal to the physical mean free path. The total particle number 
in the box is 



Pp^' 



pp(Sx)^ ' 

This gives a mean free path for the "grid point particles" of 



(C.4) 



A^ = ^='-^ = '-lSx = A. 
No- pp(Sx^) Pp 



(C.5) 



This shows how the superparticle Monte Carlo method smoothly 
transforms to the inflated particle method when reducing the 
number of superparticles and increasing their mass. At the point 
when the superparticle fills up its grid cell, the collision proba- 
bility approaches unity inside the cell and the mean free path of 
the grid cell particles is equal to the mean free path of the physi- 
cal particles. At the same time one must only allow approaching 
particles to collide, to avoid multiple collisions inside the grid 
cell. Of course, the collision detection algorithm for these cubic 
particles is rather crude, but the geometric eff'ect of considering 
cubic rather than spherical particles is minor. 



